import pylab as pl #导入绘图模块
from matplotlib.pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']   #显示中文
mpl.rcParams['axes.unicode_minus']=False       #显示负号
import pandas as pd
import numpy as np
from PyEMD import EMD, Visualisation
from dtaidistance import dtw

from scipy.signal import argrelextrema

#进行样条差值
import scipy.interpolate as spi
import matplotlib.pyplot as plt

import math
import numpy as np
import matplotlib.pyplot as plt
def moving_average(interval, windowsize):
    window = np.ones(int(windowsize)) / float(windowsize)
    re = np.convolve(interval, window, 'same')
    # re = np.convolve(interval, window, 'full')
    # re = np.convolve(interval, window, 'valid')
    return re
data = pd.read_excel('H033_L1_QDS_DCE_20230823153615675.xlsx')
coilid = 136999573300
# coilid = 137049470100
data1 = data[data['coilid'] == coilid]
data1['temp1_2'] = data1['temp1_2'].astype(float)
data11 = data1[data1['temp1_2'] >= 400]
data11 = data11.reset_index(drop=True)
df0 = data11.copy()
df0['time'] = df0['pos'].astype(float)
df0['temp1_2'] = df0['temp1_2'].astype(float)
# df0 = df0.head(100)
timestamp_array = df0['time'].values
value_array = df0['temp1_2'].values
t = timestamp_array
y = value_array

N = 5
#前5本身后4取平均
y_av = moving_average(y, N)
#解决边界点
if N % 2 == 0:
    for i in range(1, int(N / 2) + 1):
        print(i)
        y_av[i-1] = y_av[i-1] * N / (int(N / 2) + i - 1)
    for j in range(1, int(N / 2)):
        print(j)
        y_av[-j] = y_av[-j] * N / (int(N / 2) + j)
    print("偶数")
else:
    for i in range(1, int(N / 2) + 1):
        print(i)
        y_av[i-1] = y_av[i-1] * N / (int(N / 2) + i)
    for j in range(1, int(N / 2) + 1):
        print(j)
        y_av[-j] = y_av[-j] * N / (int(N / 2) + j)
    print("奇数")
df0['new_value'] = y_av



x = [float(v) for v in df0['new_value']]
lx = len(x)
print("共有%d点数据" % (len(x)))
sampling_rate = len(x)  # 取样频率(来自传感器说明书)
fft_size = len(x)  # FFT处理的数据样本数
x = np.asarray(x)
data111 = df0[515:665]
# data111 = df0[285:385]

data111['value'] = data111['new_value'].astype(float)
x1 = data111.value.tolist()
l = len(x1)
df_out = pd.DataFrame(columns=['i', 'dtw'])
dict = {}
print(l)



# coilid = 136999573300
coilid = 137049470100

data1 = data[data['coilid'] == coilid]
data1['temp1_2'] = data1['temp1_2'].astype(float)
data11 = data1[data1['temp1_2'] >= 400]
data11 = data11.reset_index(drop=True)
df0 = data11.copy()
df0['time'] = df0['pos'].astype(float)
df0['temp1_2'] = df0['temp1_2'].astype(float)
# df0 = df0.head(100)
timestamp_array = df0['time'].values
value_array = df0['temp1_2'].values
t = timestamp_array
y = value_array

N = 5
#前5本身后4取平均
y_av = moving_average(y, N)
#解决边界点
if N % 2 == 0:
    for i in range(1, int(N / 2) + 1):
        print(i)
        y_av[i-1] = y_av[i-1] * N / (int(N / 2) + i - 1)
    for j in range(1, int(N / 2)):
        print(j)
        y_av[-j] = y_av[-j] * N / (int(N / 2) + j)
    print("偶数")
else:
    for i in range(1, int(N / 2) + 1):
        print(i)
        y_av[i-1] = y_av[i-1] * N / (int(N / 2) + i)
    for j in range(1, int(N / 2) + 1):
        print(j)
        y_av[-j] = y_av[-j] * N / (int(N / 2) + j)
    print("奇数")
df0['new_value'] = y_av



x = [float(v) for v in df0['new_value']]
lx = len(x)
print("共有%d点数据" % (len(x)))



#
#
#
#
#
#
#
# def __cal_red_blue(x):
#     if x.pos>=515 and x.pos<515+l:
#         rst = 'black'
#     elif x.pos>=663 and x.pos<663+l:
#         rst = 'red1'
#     else:
#         rst = 'blue'
#
#
#
#     # if x.pos>=285 and x.pos<285+l:
#     #     rst = 'black'
#     # elif x.pos>=111 and x.pos<111+l:
#     #     rst = 'red1'
#     # elif x.pos>=383 and x.pos<383+l:
#     #     rst = 'red2'
#     # elif x.pos>=478 and x.pos<478+l:
#     #     rst = 'red3'
#     # elif x.pos>=600 and x.pos<600+l:
#     #     rst = 'red4'
#     # else:
#     #     rst = 'blue'
#     return rst
# df0['mark'] = df0.apply(lambda x: __cal_red_blue(x), axis=1)
#
# data1111 = df0[df0['mark']=='red1']
# data9999 = df0[df0['mark']=='black']
#
#
#
# # data1111 = df0[df0['mark']=='red1']
# # data2222 = df0[df0['mark']=='red2']
# # data3333 = df0[df0['mark']=='red3']
# # data4444 = df0[df0['mark']=='red4']
# # data9999 = df0[df0['mark']=='black']
# x00 = df0['pos'].values
# y00 = df0['new_value'].values
#
# x11 = data1111['pos'].values
# y11 = data1111['new_value'].values
# x99 = data9999['pos'].values
# y99 = data9999['new_value'].values
# # x11 = data1111['pos'].values
# # y11 = data1111['new_value'].values
# # x22 = data2222['pos'].values
# # y22 = data2222['new_value'].values
# # x33 = data3333['pos'].values
# # y33 = data3333['new_value'].values
# # x44 = data4444['pos'].values
# # y44 = data4444['new_value'].values
# # x99 = data9999['pos'].values
# # y99 = data9999['new_value'].values
# # plt.cla()
# plt.figure(figsize=(20, 8))
# plt.plot(x00, y00, linewidth =5.0, color = 'b')
# plt.plot(x11, y11, linewidth =5.0, color = 'r')
# plt.plot(x99, y99, linewidth =5.0, color = 'black')
# # plt.plot(x00, y00, linewidth =2.0, color = 'b')
# # plt.plot(x11, y11, linewidth =2.0, color = 'r')
# # plt.plot(x22, y22, linewidth =2.0, color = 'r')
# # plt.plot(x33, y33, linewidth =2.0, color = 'r')
# # plt.plot(x44, y44, linewidth =2.0, color = 'r')
# # plt.plot(x99, y99, linewidth =2.0, color = 'black')
# plt.ylim(550, 650)
# # plt.ylim(550, 600)
# plt.xlabel('Pos')
# plt.ylabel('Value')
# plt.show()
#
# print('finish')





for i in range(0,lx-l):
    print(i)
    data222 = df0[i:i+l]
    data222['value'] = data222['new_value'].astype(float)
    x2 = data222.value.tolist()
    distance = dtw.distance(x1, x2)
    dict['i'] = i
    dict['dtw'] = distance
    new_row = pd.Series(dict)
    df_out = df_out.append(new_row, ignore_index=True)
    x11 = data222['pos'].values
    y11 = data222['temp1_2'].values

    x22 = data222['pos'].values
    y22 = data222['value'].values

    plt.cla()
    plt.figure(figsize=(14, 8))
    plt.plot(x11, y11, 'b')
    plt.plot(x22, y22, 'r')
    # plt.ylim(550, 650)
    plt.ylim(550, 600)
    plt.xlabel('Pos')
    plt.ylabel('Value')
    # png_name = 'D:/repos/sicost/picture_09012/'+str(i)+'.png'
    png_name = 'D:/repos/sicost/picture_09013/'+str(i)+'.png'
    plt.savefig(png_name)



# writer = pd.ExcelWriter('dtw_09012.xlsx')
writer = pd.ExcelWriter('dtw_09013.xlsx')
df_out.to_excel(writer, sheet_name='Sheet1', index=False)
writer.save()




print(data.iloc[0,1])
print(type(data.iloc[0,1]))
def get_picture(coilid):
    data1 = data[data['coilid'] == coilid]
    data1['temp1_2'] = data1['temp1_2'].astype(float)
    data11 = data1[data1['temp1_2'] >= 400]
    data11 = data11.reset_index(drop=True)

    x = [float(v) for v in data11['temp1_2']]
    print("共有%d点数据" % (len(x)))
    sampling_rate = len(x)  # 取样频率(来自传感器说明书)
    fft_size = len(x)  # FFT处理的数据样本数
    x = np.asarray(x)
    xt = x - np.mean(x)  # 去直流分量

    tstemp = np.linspace(0, int(1e6 / sampling_rate * fft_size), fft_size) / 1e3
    # # pl.cla()
    # pl.figure(figsize=(14, 8))
    # title = str(coilid) + '时间趋势图'
    # pl.title(title)
    # pl.ylabel("温度")
    # pl.xlabel("时间(单位ms)")
    # # pl.plot(tstemp, xt)
    # pl.plot(tstemp, x)
    # # pl.savefig(title + '.png')
    # # pl.show()
    # list_x = xt.tolist()
    # print('the type of list_x is %s' % (type(list_x)))
    # # 下一行注释掉，横坐标为时间ms，否则为数字标号
    # # tstemp = np.arange(len(xt))

    from scipy.fftpack import fft, fftfreq
    fft_series = fft(x)
    # fft_series = fft(data["value"].values)
    power = np.abs(fft_series)
    sample_freq = fftfreq(fft_series.size)

    pos_mask = np.where(sample_freq > 0)
    freqs = sample_freq[pos_mask]
    powers = power[pos_mask]

    top_k_seasons = 3
    # top K=3 index
    top_k_idxs = np.argpartition(powers, -top_k_seasons)[-top_k_seasons:]
    top_k_power = powers[top_k_idxs]
    fft_periods = (1 / freqs[top_k_idxs]).astype(int)

    print(f"top_k_power: {top_k_power}")
    print(f"fft_periods: {fft_periods}")
    print('finish')

    from statsmodels.tsa.stattools import acf

    # Expected time period
    for lag in fft_periods:
        # lag = fft_periods[np.abs(fft_periods - time_lag).argmin()]
        acf_score = acf(x, nlags=lag)[-1]
        print(f"lag: {lag} fft acf: {acf_score}")
    print('finish')

    # expected_lags = np.array([timedelta(hours=12) / timedelta(minutes=5), timedelta(days=1) / timedelta(minutes=5),
    #                           timedelta(days=7) / timedelta(minutes=5)]).astype(int)
    # for lag in expected_lags:
    #     acf_score = acf(data["value"].values, nlags=lag, fft=False)[-1]
    #     print(f"lag: {lag} expected acf: {acf_score}")





#####上下包络线
    # datax = x
    # index = list(range(len(datax)))
    # # 获取极值点
    # d = np.diff(datax)
    # d1, d2 = d[:-1], d[1:]
    # indmin = np.nonzero(np.r_[d1 * d2 < 0] & np.r_[d1 < 0])[0] + 1
    # indmax = np.nonzero(np.r_[d1 * d2 < 0] & np.r_[d1 > 0])[0] + 1
    # # When two or more points have the same value
    # if np.any(d == 0):
    #     imax, imin = [], []
    #     bad = d == 0
    #     dd = np.diff(np.append(np.append(0, bad), 0))
    #     debs = np.nonzero(dd == 1)[0]
    #     fins = np.nonzero(dd == -1)[0]
    #     if debs[0] == 1:
    #         if len(debs) > 1:
    #             debs, fins = debs[1:], fins[1:]
    #         else:
    #             debs, fins = [], []
    #     if len(debs) > 0:
    #         if fins[-1] == len(datax) - 1:
    #             if len(debs) > 1:
    #                 debs, fins = debs[:-1], fins[:-1]
    #             else:
    #                 debs, fins = [], []
    #     lc = len(debs)
    #     if lc > 0:
    #         for k in range(lc):
    #             if d[debs[k] - 1] > 0:
    #                 if d[fins[k]] < 0:
    #                     imax.append(np.round((fins[k] + debs[k]) / 2.0))
    #             else:
    #                 if d[fins[k]] > 0:
    #                     imin.append(np.round((fins[k] + debs[k]) / 2.0))
    #     if len(imax) > 0:
    #         indmax = indmax.tolist()
    #         for x in imax:
    #             indmax.append(int(x))
    #         indmax.sort()
    #     if len(imin) > 0:
    #         indmin = indmin.tolist()
    #         for x in imin:
    #             indmin.append(int(x))
    #         indmin.sort()
    # max_peaks = indmax
    # min_peaks = indmin
    # # # 将极值点拟合为曲线
    # ipo3_max = spi.splrep(max_peaks, datax[max_peaks], k=3)  # 样本点导入，生成参数
    # iy3_max = spi.splev(index, ipo3_max)  # 根据观测点和样条参数，生成插值
    # ipo3_min = spi.splrep(min_peaks, datax[min_peaks], k=3)  # 样本点导入，生成参数
    # iy3_min = spi.splev(index, ipo3_min)  # 根据观测点和样条参数，生成插值
    # # 计算平均包络线
    # # iy3_mean = (iy3_max + iy3_min) / 2
    # # 绘制图像
    # plt.figure(figsize=(14, 8))
    # plt.plot(tstemp, datax, label='Original')
    # plt.plot(tstemp, iy3_max, label='Maximun Peaks')
    # plt.plot(tstemp, iy3_min, label='Minimun Peaks')
    # plt.legend()
    # plt.xlabel('time')
    # plt.ylabel('value')
    # title = str(coilid) + '时间趋势图及上下包络线'
    # plt.savefig(title + '.png')


####EMD分解IMF
    # t = np.array(tstemp)
    # # 提取imfs和剩余信号res
    # emd = EMD()
    # emd.emd(x)
    # imfs, res = emd.get_imfs_and_residue()
    # # 绘制 IMF
    # vis = Visualisation()
    # vis.plot_imfs(imfs=imfs, residue=res, t=t, include_residue=True)
    # # # 绘制并显示所有提供的IMF的瞬时频率
    # # vis.plot_instant_freq(t, imfs=imfs)
    # vis.show()

#####频域图
    # am = np.fft.fft(xt)  # 对希尔伯特变换后的at做fft变换获得幅值
    # am = np.abs(am)  # 对幅值求绝对值（此时的绝对值很大）
    # am = am / len(am) * 2
    # am = am[0: int(len(am) / 2)]
    # freq = np.fft.fftfreq(len(xt), d=1 / sampling_rate)  # 获取fft频率，此时包括正频率和负频率
    # freq = freq[0:int(len(freq) / 2)]  # 获取正频率
    #
    # pl.figure(figsize=(14, 8))
    # title = str(coilid) + '去直流频域图'
    # pl.title(title)
    # pl.ylabel("振幅")
    # pl.xlabel("频率")
    # pl.plot(freq, am)
    # # pl.savefig(title + '.png')
    # # pl.show()
get_picture(136999573200)
get_picture(136999573300)
get_picture(136999573400)
get_picture(136999573500)
get_picture(136999573600)
get_picture(137049470100)

print('finish')


